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ABSTRACT 

The interaction of a fast wind with a spherical Asymptotic Giant Branch (AGB) wind is thought to 
be the basic mechanism for shaping Pre-Planetary Nebulae (PPN) and later Planetary Nebulae (PN). 
Due to the large speed of the fast wind, one expects extended X-ray emission from these objects, but 
X-ray emission has only been detected in a small fraction of PNs and only in one PPN. Using numerical 
simulations we investigate the constraints that can be set on the physical properties of the fast wind 
(speed, mass-flux, opening angle) in order to produce the observed X-ray emission properties of PPNs 
and PNs. We combine numerical hydrodynamical simulations including radiative cooling using the 
code FLASH with calculations of the X-ray properties of the resulting expanding hot bubble using 
the atomic database ATOMDB. In this first study, we compute X-ray fluxes and spectra using one- 
dimensional models. Comparing our results with analytical solutions, we find some agreements and 
many disagreements. In particular, we test the effect of different time histories of the fast wind on the 
X-ray emission and find that it is determined by the final stage of the time history during which the fast 
wind velocity has its largest value. The disagreements which are both qualitative and quantitative in 
nature argue for the necessity of using numerical simulations for understanding the X-ray properties of 
PNs. The X-ray luminosity in the 0.2 - 10 keV range (1.24 - 62 A) covered by the CHANDRA/^C/S' 
instrument shows a very slow decrease with time over the range of evolutionary ages explored in our 
models (up to 3750 years). Furthermore, investigating the emission in other wavelengths ranges, we 
find that most of the luminosity emerges at longer wavelengths (A > 140 A) from the cooler outer 
edge of the hot bubble. We apply our spherical models to the objects BD-|-30°3639 and NGC 40. We 
find that the model values of the X-ray temperature and luminosity for these objects are significantly 
higher than observed values and discuss several mechanisms for resolving the discrepancies. 
Subject headings: circumstellar matter — hydrodynamics — ISM: jets and outflows — planetary 
nebulae: general — stars: mass loss — X-rays: ISM 



1. INTRODUCTION 

The shaping of Pre-planetary nebulae (PPN) and Plan- 
etary Nebulae (PN) is believed to result from the inter- 
action of a fast, coUimated post- AGB wind plowing into 
the slow, dense wind emitted during the AGB phase, fol- 
lowed by an isotropic tenuous wind during the PN phase 
(Sahai & Trauger 1998). These interacting winds result 
in an expanding shell of shocked matter which forms the 
PPN and later the PN. Due to the large speed of the 
fast wind (from few x 100 km s~^ up to 2000 km s^^), 
one expects extended X-ray emission to be produced in 
PPNs (where the shaping mainly occurs) and PNs. 

However, X-ray emission was only detected in 3 of 60 
PNs observed with ROSAT, followed by a few more from 
CHANDRA and XMM (e.g. Guerrero et al. 2005). In 
the case of PPNs, there is so far only one confirmed X- 
ray detection (Sahai et al. 2003), although many have 
been observed with CHANDRA. X-ray emission should 
be present in principle in all of these systems ~ hence the 
low detection rate implies that it is obviously below the 
detection limit. 

The general problem of understanding the formation 
and shaping of PNs has been addressed analytically and 
numerically (see Balick & Frank 2002, and references 
therein). In principle, one of the most direct probes of 
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the fast wind and the interaction process which drives PN 
formation is the X-ray emission. It may help in answering 
several important questions - e.g. what are the physical 
properties of the fast wind (speed, mass-flux, opening an- 
gle - and the time histories of these parameters). Such 
studies can also help in identifying the nebular regions 
that are mainly responsible for the X-ray emission for 
comparison with recent X-ray images of PNs obtained 
with CHANDRA and XMM (e.g. Kastner et al. 2000, 
2003; Guerrero et al. 2005). 

X-ray emission as a probe was first exploited by Zhekov 
& Perinotto (1996), who derived an analytical, self- 
similar, spherically symmetric model, taking into account 
the effects of thermal conductivity. Very recently, Akashi 
et al. (2006, hereafter ASB06) presented new analytical 
models with radiative cooling, but without heat conduc- 
tion. Although analytical modeling is potentially a pow- 
erful tool for understanding of the X-ray emission from 
PNs and its dependence on fundamental physical param- 
eters, it necessitates making certain assumptions which 
may or may not be completely valid. For example, a ma- 
jor problem with such analytical models is, that they are 
not able to treat all the effects (hydrodynamics, radiative 
cooling, heat conduction) self-consistently. In this paper, 
we have therefore performed numerical simulations using 
the code FLASH in which we vary the basic parameters 
of the fast and slow wind over an extensive parameter 
grid and compute the X-ray emission as a function of 
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these parameters. A major goal of our work here is to 
test the main results from ASB06's analytical modeling 
and verify their general conclusions. We have therefore, 
as in ASB06, not included heat conduction in our treat- 
ment. However, we discuss the heat conduction process 
and how it may be inhibited by magnetic fields which 
are likely to be present in PNs. For example, the lack of 
spherical symmetry in most PNs has been explained by 
some authors as resulting from the presence of magnetic 
fields (e.g. Garcia-Segura et al. 2005). This paper, in 
whic;li we have carried out spherically symmetric simula- 
tions of interacting winds, is a first in a series of papers 
in our quest to understand X-ray emission from PPNs 
and PNs. 

The rest of the paper is organized as follows. In §2, 
we describe our numerical model. In §3, we present the 
results of our one-dimensional simulations including the 
structure and X-ray properties of the expanding bubble 
and their dependency on the parameters of the fast wind. 
In §4, new effects from the computation of a spherically 
symmetrical model on a two-dimensional grid are shown. 
In §5, we apply our models to two PNs, BD-|-30°3639 and 
NGC 40, and examine several mechanisms for resolving 
the discrepancies between the models and observed val- 
ues of the X-ray temperature Tx and luminosity for 
these objects. We compare in detail our results with 
those of ASB06 in §6 and discuss the still unsolved prob- 
lem of fitting Tx and in PNs in §7. Finally in §8, we 
present our conclusions. 

2. THE NUMERICAL METHOD 

2.1. Initial and boundary conditions 

Using FLASH, we solve the following set of the differ- 
ential equations of ideal hydrodynamics 



and 



ff + V(pv) = 



9(pv) 

dt 



V (pv (g) v) = -Vp 



de 

'dt 



+ V (e v) = -p V V - A(T) 

P=(7- 



(1) 

p is the gas density, p the pressure, e the internal energy 
density, v the velocity and 7 the ratio of the specific heats 
at constant pressure and volume which is set to 7 = 5/3. 
We include radiative losses using the cooling function 
A(r) of Sutherland & Dopita (1993) in the temperature 
range between 10^ and 10®-^ K and the low temperature 
function of Dalgarno & McCray (1972) for T < 10* K. 

A spherically symmetric AGE wind is implemented 
with 

Ms 
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and 
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The temperature of the AGB wind is set to 10 K. In this 
r~^-density profile, we let a coUimated fast wind (CFW) 
plow from the inner radial boundary with the following 
parametrization: 
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1 2^ 
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V{ = V{ exp 



1 2^ 
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9f is the opening angle of the CFW. K is a normalization 
factor calculated by 



K 



it/2 



exp 



sin 0d0. 
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For a temperature of the CFW, we assume 10** K. The 
simulations are performed in spherical coordinates with 
the dimensions of r = 5 x 10^'^ - 2.5 x 10^^ cm and 

9 = 0- 7r/2 in the case of the two-dimensional runs. 
Mirror boundaries are used along the symmetry axis and 
the equatorial plane and outflow conditions at the upper 
radial boundary. 

2.2. Model parameters 

The models which are presented in this paper are 
spherical ones without the 9 dependency in eqns. 4 - 
5 and K = 1. Preliminary results of models with 9 de- 
pendency (i.e. incorporating collimated fast winds) were 
published in (Stute & Sahai 2006). Detailed results are 
deferred to a future publication. The parameters of the 
models are then the mass outflow rate and velocity of 
the fast and the slow wind, respectively. For these pa- 
rameters we adopted the same ranges as in the models 
of ASB06 for ease of comparison (see Table 1). Thus we 
performed three sets of runs with Mg of 3 x 10^^ (runs 
A), 7 X 10"*' (nms B) and 3 x IQ-"^ Mq yr^^ (runs C) 
and within each set we varied the velocity of the fast wind 
from 300 to 700 km s'^ (e.g. the 5 in 05 means Vf = 500 
km s"""^). The velocity of the slow wind was set to 10 
km s~-^ and Mf was chosen such that M{ Vf = Mg Vg. We 
expanded these sets while applying our models to obser- 
vations. 

2.3. Calculating the X-ray properties 

Using the radial density and temperature profiles from 
the hydrodynamical simulations, we determine the ex- 
pected X-ray flux. We have used the atomic database 
ATOMDB with IDL including the Astrophysical Plasma 
Emission Database (APED) and the spectral models 
output from the Astrophysical Plasma Emission Code 
{APEC) (Smith et al. 2001) to calculate the emissivity. 
The default abundances in ATOMDB, i.e. 14 elements 
(H, He, C, N, O, Ne, Mg, Al, Si, S, Ar, Ca, Fe, Ni) with 
solar abundances of Anders & Grevesse (1989) are used. 
The energy range is divided into bins of 0.01 keV. We 
compute the spectrum and the total flux in X-rays as a 
function of evolutionary time for each of our models. 

X-ray emission in an interacting winds model, in prin- 
ciple, will arise from gas with a large range of tempera- 
tures above lO"' K. We calculate the X-ray emission be- 
tween 0.01 - 10 keV. The lower energy boundary is deter- 
mined by the limitations of ATOMDB. We have divided 
this range into three bins. The high energy bin (0.2 - 

10 keV or 1.24 - 62 A) represents the energy range cov- 
ered by the A (7/5 instrument and (almost) that of HETG 
(0.4 - 10 keV) on CHANDRA. We define a medium en- 
ergy bin (0.09 0.2 keV or 62 137 A), since the energy 
ranges of the EPIC instrument on XMM (0.1 - 12 keV) 
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TABLE 1 



run 


Mf 


V{ 


Ms 




Vci 






J 


fASB 


^cd.sim 


A3 


1 X 10"^ 


300 


3 X 10-"^ 


10 


21.54 


27.51 


9.58 


0.033 


28.6 


25.16 


A5 


6 X 10"** 


500 


3 X 10-'^ 


10 


25.54 


32.62 


9.8 


0.02 


32.6 


28.47 


A7 


4.33 X 10"** 


700 


3 X 10"^ 


10 


28.67 


36.62 


9.66 


0.014 


35.6 


30.92 



B3 


2.33 X 10" ' 


300 


7 > 


< 10" 


6 


10 


21.53 


27.51 


9.58 


0.033 


28.6 


23.97 


B4 


1.75 X 10"^ 


400 


7 X 10" 


6 


10 


23.71 


30.28 


9.75 


0.025 


30.8 


25.70 


B5 


1.4 X lO"'' 


500 


7 X 10- 


6 


10 


25.-54 


32.62 


9.8 


0.02 


32.6 


26.96 


B6 


1.17 X 10"^ 


600 


7 > 


< 10- 


6 


10 


27.17 


34.7 


10.03 


0.017 


34.2 


28.10 


B7 


1 X 10"^ 


700 


7 > 


< 10- 


6 


10 


28.58 


36.49 


9.66 


0.014 


35.6 


29.10 


C3 


1 X 10"" 


300 


3 > 


< 10- 


5 


10 


21.54 


27.51 


9.58 


0.033 


28.6 


22.60 


C5 


6 X 10"^ 


500 


3 > 


< 10- 


5 


10 


25.54 


32.62 


9.8 


0.02 


32.6 


25.05 


C7 


4.28 X 10"^ 


700 


3 > 


< 10" 


.5 


10 


28.57 


36.47 


9.66 


0.014 


35.6 


26.88 


AlO 


3 X 10"" 


1000 


3 X 10" 


i; 


10 


32.18 


41.10 


9.90 


0.01 




36.56 


BIO 


7 X 10"* 


1000 


7 X 10- 


6 


10 


32.18 


41.10 


9.90 


0.01 




33.43 


D7 


1 X 10"'^ 


700 


7 X 10- 


5 


10 


28.58 


36.49 


9.66 


0.014 




25.99 


ElO 


2.4 X 10-6 


1000 


5 X 10- 


4 


10 


25.19 


32.17 


6.88 


0.007 




27.03 


B5 (2D) 


1.4 X IQ- ' 


500 


7 X 10- 


6 


10 


25.54 


32.62 


9.8 


0.02 




22.4 



Note. — The first five columns give the parameters of our models: the name of the run, the 
mass outflow rate and velocity of the fast and the slow wind, respectively. The next four columns 
give results of the analytical solution of Koo & McKee (1992a,b): the critical velocity, the velocity 
of the contact discontinuity in the energy-driven case and in the momentum-driven case and the 
constant X. The velocity used by ASB06 following a formula in Volk & Kwok (1985) is given in 
the column before last, the last column shows the measured velocity of the contact discontinuity 
iicd.sim in the simulations. 



and of LETG on CHANDRA (0.09 - 3 keV) extend to 
lower energies, and a low energy bin between 0.01 - 0.09 
keV (137 - 1180 A) which is sensitive to gas with tem- 
peratures < 10^ K. In Fig. 1, wc show how the X-ray 
emissivity varies as a function of temperature in the three 
energy bins we have defined above. It is interesting to 
note that gas with a temperature lower than 2 x 10^ K 
contributes only marginally in the high energy bin which 
is probed by the ACIS instrument on CHANDRA (Fig. 
1, bottom). 

3. ONE-DIMENSIONAL SIMULATIONS 
3.1. Structure 

The profiles of density and temperature as a function of 
radius (Fig. 2; runs B3, B5 and B7 at an age of 390 years) 
show four different regions: (i) the unshocked fast wind, 
(ii) the shocked fast wind (hot bubble), (iii) the shocked 
slow wind forming a dense shell and (iv) the unshocked 
slow wind. The reverse shock separates region (i) from 
(ii). The large jump in density and temperature between 
regions (ii) and (iii) represents the contact discontinuity. 
The outer shock separates region (iii) from (iv). Due to 
the interaction of the fast and the slow wind, both regions 
of shocked matter are heated to high temperatures, but 
the dense shell is highly radiative and cools quickly. The 
evolution of the shell normally has the following stages: 
first there is free expansion, if cooling sets in there is then 
the radiative, momentum-driven bubble and finally the 
adiabatic, energy-driven bubble. 

Depending on the velocity of the fast wind, namely if 
it is larger or smaller than a critical velocity 

V67rpoi/ \3 Ms J 

the end stage is then the momentum-driven or energy- 
driven bubble (Koo & McKee 1992a,b), where £ = 
l/2M{Vf is the kinetic luminosity of the fast wind. This 
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Fig. 1. — Top: X-ray emissivity between 0.01 - 10 keV (solid), 
0.09 - 10 keV (dash-dotted) and 0.2 - 10 keV (dashed) derived from 
the ATOMDB as a function of the plasma temperature; bottom: 
X-ray emissivity in the three different energy bins 0.01 — 0.09 keV 
(solid), 0.09 - 0.2 keV (dash-dotted) and 0.2 - 10 keV (dashed) 

critical velocity is independent of the cooling function 
in the case of a constant mass outflow rate of the fast 
wind and a r"^ density profile (Koo & McKee 1992b). 
As the fast wind is always faster than in our mod- 
els, the bubble is expected to go from free expansion 
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run B7 
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run 87 ■ 
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Fig. 2. — Plots of the density and temperature at an age of 390 
years in the one-dimensional runs B3 (dotted), B5 (solid) and B7 
(dashed), i.e., Ms = 7 X IQ-'^ M© yr^^ and V{ is 300, 500 and 700 
km s~^, respectively; the position of the dense shell (ds) and hot 
bubble (hb) arc indicated for model B5 in the density plot 

phase directly to the energy-driven phase. This transi- 
tion occurs in our models within a short period of time 

(about few years), after which the velocity of the contact 
discontinuity reaches a stable value characteristic of an 
energy-driven shell. In the momentum-driven case, the 
position of the shell would be (Koo & McKee 1992a, b) 



R: 



M . 
■cd ■ 



I 



V{ t 



1/2 



with 



L 



2 TT poi 



1/2 



1/2 



(8) 



(9) 



The resulting velocity of about 10 km s^^ (Table 1) are 
significantly lower than those found in our simulations. 
In the energy-driven, adiabatic case, the analytical solu- 
tions of the position of the shell, i?cd> is 



pE . 



■ rad 



1/3 



t 



. 3poi 

' 2 7rrrad^Mft;f^t;s 
3M« 



1/3 



(10) 



with 47rrrad^ = 4.165 (Koo & McKee 1992a,b). Here 
the ambient shock is assumed to be adiabatic. Assum- 
ing a radiative ambient shock, 4 tt Fj-ad C = 2.0 and the 

velocity of the shell reduces to v^r- is the fraction 

of the injected energy in the bubble which is radiated 



away, ^ is a numerical constant of order unity (Koo & 
McKee 1992b). The velocity of the contact discontinuity 
Wcdjsim in our simulations is within the range of velocities 
between the limiting values given by the radiative and 
adiabatic approximations - Vcr and Rf^/t - for the dense 
sheh of Koo & McKee (1992a,b) (see Table 1). 

The density and temperature of the hot btibble are in 
good agreement with results from the self-similar formu- 
lation of Chevalier & Imamura (1983), although radiative 
cooling emphasizes the density and temperature jumps 
close to the contact discontinuity. For example in our 
model B5, the density is about 3 x 10^^^ g cm^^, i.e. a 
number density of about 200 cm^'', and the temperature 
is 5 X 10^ K. Our density is slightly lower than the value 
derived from the self-similar formulation of 350 cm~^ - 
this leads to a reduced cooling efficiency, and thus our 
temperature is slightly higher than the value from the 
self-similar model of 3.4 x 10^ K. 

The temperature of the hot bubble is given in Table 2. 
We calculated an X-ray weighted average temperature as 
defined in ASB06. It is only a function of the velocity 
of the fast wind, ranging from 1.5 x 10^ K for Vf — 300 
km to 2.2 X 10^ K for Vf = 1000 km s'^ The mass 
loss rate of the slow wind is fairly unimportant, i.e., the 
temperature is almost identical, e.g., in the models A3, 
B3 and C3. The temperature is roughly given by 



Thb = 1.4 X lO^'K 



300 kms- 



(11) 



in which the dependence on only Vf is as expected from 
the Rankine-Hugoniot conditions (e.g. Landau & Lifshitz 
1959). As the bubble evolves, the temperature of the hot 
bubble remains fairly constant with time. 

3.2. X-ray properties 

We now present the properties of the X-ray emission 
from our models and their dependence on the physical 
parameters of the fast wind. 

3.2.1. The X-ray luminosities 

Several clear trends in the X-ray luminosity are visible 
between the different sets of simulations. We show the 
X-ray luminosity in the high energy bin between 0.2 - 10 
keV, ix.ACiSj as a function of the evolution time for all 
models in Fig. 3. Comparing models from different sets 
with fixed Mg which have the same fast wind velocity, 
we find that the higher the mass outflow rate of the fast 
wind is, the higher Lx.Acis is. For example, at an age of 
470 years, in runs A3, B3 and C3 Lx,ACiS is 0.1 x 10^^ erg 
s~\ 0.4x 10^^ erg s~^ and 2.4x 10^^ erg s~\ respectively. 
In these runs, the fast wind velocity is 300 km s"""^ and 
Mf is 10"'^, 2.33 X 10"'^ and 10-*^ Mq yr~\ respectively. 

Comparing the curves within each panel (each panel 

represents a fixed Ms and thus a constant Mf V{), ix.ACis 
decreases with increasing fast wind velocity. This is be- 
cause the correspondingly lower mass outflow rate of the 
fast wind^ leads to lower densities in the shocked fast 
wind and therefore to reduced X-ray emission. 

We can also investigate the dependence of ix.AciS on 
Vf, when Mf is constant. Using pairs of our runs, as 

1 A result of the assumption Mf tif = Ms f s 
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TABLE 2 



run 


Red ( 10i« cm ) 


-Rcd/^ASB ( yrs ) 


Lasb ( erg s 1 ) 


i-x.ACIS ( 


erg s 1 ) 


Thb ( K ) 


Ad 


3.4 


377 


i.O X iU 


4.33 X 10^1 


1.5 X 10 


A5 


3.9 


379 


2. X lU 


2.48 X 10^1 


4.5 X lU 


A7 


4.3 


383 


1.3 X 10 


1.20 ; 


< 103i 


8.8 X 10 


B3 


3.2 


355 


6.9 X 10^^ 


1.64 ; 


< 10^^ 


1.5 X 10" 


B4 


3.5 


360 




1.48 ; 


< 1032 


3.0 X 10*^ 


B5 


3.7 


361 


1.1 X 10^2 


9.48 X 10^1 


4.9 X 10*^ 


B6 


3.8 


352 




6.84 X 10-" 


7.1 X 10*5 


B7 


4.0 


356 


7.6 X 10-" 


4.31 ; 


< 1031 


9.6 X 10^' 


C3 


3.0 


332 


< 2. X lO^i 


1.01 ; 


< 1033 


1.7 X 10" 


C5 


3.4 


330 


1.4 X 10^3 


5.03 ; 


< 1032 


5.3 X 10" 


C7 


3.7 


329 


1.3 X 10^3 


2.40 ; 


< 1032 


1.1 X 10^ 


AlO 


4.7 






2.84 ; 


< 103" 


1.8 X 10' 


BIO 


4.4 






1.00 X 10-" 


2.0 X 10^ 


D7 


4.0 






9.05 X 10^2 


1.1 X 10^ 


ElO 


3.0 






2.34 X 10^3 


2.2 X 10'^ 



Note. — The X-ray luminosity Lx.ACiS in ' 
years in our simulations and results Lasb from 
size in our models and those of ASB06, since 1 
different in our models compared to ASB06); X-: 
Tjjb as defined in ASB06 at an age of 470 years 

e.g. runs A3 and B7 (in which Mf = lO""^ Mq yr"!) 
or C3 and D7 (Mf = 10"^ M© yr"!), we can see that 
the luminosities arc almost the same in each pair (Ta- 
ble 2), even though, due to the different velocities, the 
temperatures in the hot bubble are very different. This 
is a somewhat surprising result, but which can be un- 
derstood in terms of the important physical quantities 
determining the X-ray luminosity. These are the vol- 
ume, density and temperature of the hot bubble. The 
volumes of the hot bubble in run C3 and D7 are about 
6.84 X lO''^ cm3 and 2.2 x 10^° cm^, respectively (derived 
from Red = 3.0 x 10^^ cm and Rrs = 2.2 x 10^^ cm for C3 
and Red = 3.9 x 10^^ cm and R,s = 2.2 x 10^^ cm for D7). 
The mean (radially averaged) density in the hot bubble 
for C3 (D7) is 1.66 x IQ-^i g cm'^ (6.43 x lO'^^ g cm'^), 
the emissivity at the mean temperature of the hot bubble 
- 1.73 X 10^ K (1.06 X 10^ K) - is 1.44 x IQ-^^ erg s"^ 
cm^ (2.63 X 10^^''* erg s^^ cm''), thus the luminosity is 
9.7 x 10^2 erg s^^ (8.6 x IQ-''^ s-i)_ These estimates 
are in good agreement with the more accurate values of 
the luminosity in Table 2. In the other pair A3 and B7, 
the differences between the estimates using radially aver- 
aged quantities and the values in Table 2 are somewhat 
larger the luminosities in nms A3 and B7 are 4.9 x 10^^ 
erg s^^ and 3.9 x 10'^^ erg s^^, respectively, compared to 
4.3 X 10^^ erg s~^ in both cases in Table 2. This is be- 
cause the density and temperature gradients across the 
hot bubble are larger than in the first pair and therefore 
also the errors, we introduce by the averaging process. 

In summary, the weak dependence of the X-ray lumi- 
nosity, Lx,ACiSi on Vf has two reasons: i) the emissivity 
in the high energy bin is fairly flat over the range of tem- 
peratures in our models (Fig. 1) and ii) the volume of 

the hot bubble depends weakly on V[ {Red ~ ^'^) due to 
our assumption of Mf Vf = Ms Vg- In Fig. 4, we illustrate 
that ix.ACiS depends much more strongly on Mf than on 
Vf by plotting the values for all computed models (there- 
fore the full range of Vf). The data can be reasonably 

well fitted by the curve ix^Acis ^ Mj- the scatter of 
the individual data points above and below this curve lies 



high-energy bin 0.2 - 10 keV at an age of 470 
B06 (from epochs where the bubbles have equal 
velocity of the contact discontinuity is slightly 
weighted average temperature of the hot bubble 

within a factor of 2. Interestingly, even model ElO where 

Mf vi = O.SMg Vs fits the Lx,ACiS ~ ' relationship. 

We now show that this variation of ix,Acis with Mf 
can be partially understood by considering the dependen- 
cies of the volume {V), temperature (Thb) and density of 
the hot bubble (phb) on Mf and V{ across and within the 
different sets of models. We utilize analytical estimates 
of these dependencies based on the self-similar solutions 
without radiative cooling (see §3.1), therefore we expect 
some differences between the dependencies derived be- 
low and those actually found from our modeling (which 
takes radiative cooling effects into account). The volume 
of the hot bubble depends on the radii of the reverse 

shock (i?rs = [3Mf Wf Ws/(4M)]^/^i = ^/3/^Vst) and 
the contact discontinuity (given by eqn. 7). In the range 
of Thb in our models (1.5 x 10^ - 2 x lO'^ K), the emis- 
sivity changes only within a factor of 1.5 (Fig. 1) and 
therefore we assume it to be constant. The density of the 
hot bubble can be approximated as constant with radius 
(Fig. 2) and equal to the density of the fast wind given 
by eqn. (4) at the reverse shock scaled by a constant 
compression factor of 4, valid for an adiabatic shock. 

In general, ix.ACis is a function of the 4 parameters 
in our models - V{, Mf, Vg and Mg. Since Vg is kept 
constant in all our models and we assume Mf vt — Mg Vg, 
the number of independent parameters is only two - we 
select Vf and Mg- By holding each of these constant in 
turn one can get insights into the behavior of ix,ACis as 
a function of Mf. 

First we keep Vf constant, i.e. we compare models from 
different sets with the same V{. Then the radius of the 
contact discontinuity is constant, thus also the volume. 
Hence phb ~ M and p^^ V ~ Mf. Therefore we expect 
Lx,ACiS to vary as M^. In our models, the runs with 
the same Vf follow ix.Acis ^ M^'^^, which is primarily 
caused by a departure from the analytical estimate of the 
model (radially-averaged) density dependency, < p >~ 
MO-75. 

Next we keep Mg constant, i.e. we compare mod- 
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Fig. 3. — X-ray luminosity ix.ACiS in the high energy bin be- 
tween 0.2 - 10 keV as a function of the evolution time in the one- 
dimensional runs A3-A7 (top), B3-B7 (middle) and C3-C7 (bot- 
tom) 

els within each set, thus V{ M{ is constant. Hence 
-Red ~ ~ Mf"^/^ and therefore V ^ M^r^. The 

density phh ^ and therefore p^b^ ~ ^f' imply- 
ing Lx.Acis ^ Mf. This expectation is supported by 
the steep increase in Lx.ACiS with Mf which we see in 
sets A and B for models with the lower values of Mf, 
i.e. with the larger values of V{ (500 to 1000 km s~^). 
For our models with lower values of Wf, however, the 
increase is not as steep, due to two effects: (i) devia- 
tions from the analytical dependencies as well as (ii) the 
replacement of the radial integral for £x.acis with the 
product of radial-averages. For example, as a result of 
effect (i) in our models B3 and B4 we find R,-s ~ Mf-^^, 

Red ~ M^"■^^ V - ~ Furthermore 

the cmissivity in run B3 in 30 % lower compared to that 
in B4, equivalent to A ~ Mf°'^. The net result is that 
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Fig. 4. — X-ray luminosity ix.ACiS ™ the high energy bin be- 
tween 0.2 - 10 keV as a function of the mass loss rate of the fast 
wind Mf for all computed models at an age of 470 years (plus sign: 
V{ = 300 km s~^, asterisk: 400 km s~^, cross: 500 km s~^, dia- 
mond: 600 km s~^, triangle: 700 km s~^, square: 1000 km s~^); 

■ 3/2 

the dotted line represents the curve Lx,AClS ^ ■ Iii s-H models 

except for run ElO (the model with the highest Mf) we assumed 
Mf Vf = Ms v^. The solid lines join models within each set of con- 
stant Ms . The relation appears to hold roughly for ages of PNs up 
to ~2000 years (unjoined plus signs: runs A3, B3, and C3 at an 
age of 1900 years). 

-^x,ACis p'i^.V K M}-^ . As a result of effect (ii), the 
radial integral becames increasingly overestimated by the 

product of radially-averaged quantities for models with 
the lowest values of Vf, where the radial gradients of these 
quantities are the largest (e.g. compare run B3 and B5 
in Fig. 2). Hence the expected dependence of Lx,ACiS on 
Mf will be even shallower as observed. 

The X-ray luminosities in our different models do not 
change strongly with time as the bubble evolves, after 
an initial phase of increasing luminosity (e.g. a period of 
~ 100 years in run C3). The length of the initial phase of 
rising luminosity is primarily given by the cooling time 
in the hot bubble. Once a gas parcels has cooled below 
10^ K, it does not contribute to the Imninosity Lx.AciS 
of the hot bubble anymore. Within each set, the lower 
the value of Vi, the lower the temperature and the higher 
the density in the hot bubble, thus the shorter the cool- 
ing time. Hence in these case, the length of the initial 
phase is the highest. We discuss the weak dependence of 
ix.ACiS on time and the initial phase of rising luminos- 
ity and compare our results with those of ASB06 in §6. 

'3/2 

The relation ix.Acis ^ Mj. shown in Fig. 4 appears 
to hold for ages of PNs up to ^2000 years, presumably 
because of the similar weak decay of ix.ACiS with time 
in all models. 

Up to now, we only took into account luminosities in 
the high energy bin 0.2 10 keV, which can be observed 
with the ACIS instrument. An inspection of the X-ray 
luminosity emitted in the lower energy bins in our models 
may be of potential use in identifying new observational 
probes of cooler gas in the hot bubble near its interface 
with the dense shell. Plotting the X-ray cmissivity and 
luminosity in our three energy bins (e.g. for run B5), we 
find that most of the emission from the main body of the 
hot bubble, i.e. excluding its edges, is radiated in the 
high energy bin. However, the total energy radiated by 
the main body of the hot bubble in these three energy 
bins Lhb.m (1-4 x 10"^^ erg s^^) is significantly smaller 
than the rate at which mechanical energy is pumped 
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Fig. 5. — Top: Plot of the X-ray emissivity in the three energy 
bins 0.01 - 0.09 keV (solid), 0.09 - 0.2 keV (dotted) and 0.2 - 
10 keV (dashed) as a function of radius at an age of 390 years 
of the one-dimensional run B5; bottom: Cumulative plot of the 
luminosity inside a given radius as a function of this radius 

into the system by the fast wind, Lfw = 1.1 x lO'^^ erg 
s~^. This mechanical energy is partly converted into the 
mechanical energy of the dense shell, which increases 
at a rate given by L^a = 1/2 MsV^^/vs « 4.6 x 10^^ 
erg s^^, and the remaining energy goes into the hot 
bubble. The energy input into the hot bubble Lhb is 
about ifw — ids = 6 X 10^^ erg s~^, from of which 
Lad = ^/^ MsV^^/vs w 3 X lO''"^ erg is required for 
the adiabatic expansion of the hot bubble. Therefore the 
energy radiated by the hot bubble, Lhb — -f'ad = 3 x 10^^ 
erg s~^, is significantly larger than ihb,m and this differ- 
ence must be radiated by the edges. This is qualitatively 
demonstrated by a jump in the cumulative luminosity 
function, C^{r) (i.e., the luminosity emitted by material 
inside a given radius r), across the contact discontinu- 
ity. There are spikes in the emissivity in the low en- 
ergy bin at both edges of the hot bubble (Fig. 5, top) 
due to the presence of gas with intermediate tempera- 
tures < lO"'' ''' K. However, only the outer edge produces 
a jump in C'^{r) for the low energy bin, because it has a 
significantly higher emission measure (see Fig. 2, top). 

Hence the maximum fractional contribution to the X- 
ray luminosity emitted in the full energy range (0.01 - 
10 keV) lies in the low energy bin and comes, not from 
the main body of the hot (5 x 10^ K) bubble, but from 
its somewhat cooler (10^ - 10^ K) edges. The gas in the 
dense shell is too cool (< 10^ K) to emit any X-ray radi- 
ation. The structure of the hot bubble is similar in the 
high and medium energy bins, however, it is significantly 
limb-brightened in the low energy bin. 

A plot of the luminosity in the high and medium energy 
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Fig. 6. — X-ray luminosity as a function of the evolution time 
in the one-dimensional run B5 in the energy bins 0.09 — 0.2 keV 
(dashed) and 0.2 - 10 keV (solid); also plotted axe the mechanical 
energy injected per unit time by the fast wind (dotted) and the 
increase of mechanical energy per unit time of the dense shell (dash- 
dotted) 
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Fig. 7. — Effective area of the ACISS instrument on CHANDRA 

bins as a function of the evolution time (Fig. 6) shows 
that these are roughly constant with time. Note, that 
since most of the energy in the 0.01 - 0.09 keV range is 
emitted by only a few grid cells in our models, we do not 
have reliable quantitative estimates of the luminosity in 
the low energy bin. 

3.2.2. The X-ray spectra 

In this section, we investigate the X-ray spectra from 
our models and how they are affected by different physi- 
cal parameters. For this purpose, the model spectra have 
been convolved with the effective area of the ACIS S in- 
strument (Fig. 7). We assume a distance of 1 kpc. Our 
spectra are presented with a resolution of 0.01 keV. The 
default abundances in ATOMDB (see §2.3) are used. 

As the temperature of the hot bubble is only a func- 
tion of the fast wind velocity, the spectra from, e.g., the 
models A5, B5 and C5, look almost identical. The dif- 
ferences between each set of models are in the absolute 
values of the flux, not the relative fluxes of different lines 
in the spectra. 

Within each set, as the value of Vi changes, the spectra 
also change dramatically (Fig. 8). At moderate velocities 
of the fast wind of 300 km s~^, only a few lines with 
energies up to 1 keV are present. Increasing the fast 
wind velocity (and with it the temperature of the hot 
bubble) leads to the emergence of lines at higher energies, 
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SO that at and above velocities of 600 km s^^ one can see 
a peak of blended lines in the 0.8 - 1.1 keV range with 
an increasingly dominant contribution from iron lines. In 
the model BIO, the continuum contributes to an energy 
range from 0.3 to above 2 keV. 

4. TWO-DIMENSIONAL SPHERICAL MODELS 

In onc-dimcnsional simulations, which wc have dis- 
cussed so far, instabilities and turbulence cannot oc- 
cur. Instabilities may, however, be quite important in 
explaining microstructure in PNs, such as the observa- 
tions of globules in a few PNs (Huggins & Frank 2005). 
Therefore, we also calculated a two-dimensional spher- 
ical model^ with the parameters of run B5 to examine 
the presence of such instabilities and how they affect the 

^ i.e. a two-dimensional coordinate system and spherical initial 
conditions 



structure of the hot bubble and the dense shell and the 
resulting X-ray properties. 

The velocity of the shell is similar to the result of the 
one-dimensional simulation, but somewhat lower - the 
value is 22.4 km s"-'^ compared to about 27 km s^^. As 
the expansion in the energy-driven phase is mainly due to 
the pressure in the hot bubble, the reduction of the shell 
velocity indicates that the pressure is reduced in the two- 
dimensional run. The pressure in the one-dimensional 
run is 2.2 x 10^'' dyn and 1.47 x lO"'' dyn in the two- 
dimensional run. The latter is reduced by about one 
third, which leads to a reduction of the shell velocity of 
20 %, as Wsheii ~ \/P> i-e. a reduction from 27 to 22 
km s~^ as seen in the data. In Fig. 9, we show plots 
of the temperature and density as a function of radius 
at an age of 288 years of the one-dimensional run and 
an azimuthal average of the two-dimensional run - as 
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Fig. 9. — Plots of the density and temperature as a function of 
radius at an age of 288 years of the one-dimensional (solid) and 
two-dimensional (dotted) run B5 




Fig. 10. — Contour plots of the logarithm of density (top) and 
temperature (bottom) at an age of 288 years of the two-dimensional 
run B5 



the temperature is identical, but the density is lower, 
the pressure is lower in the latter. The reduction of the 
pressure also produces a smaller extent of the hot bubble 
(i.e., the radial distance from the reverse shock to the 
contact discontinuity) . 

The transition from a one-dimcnsional to a two- 
dimensional model introduces instabilities (Rayleigh- 
Taylor- instabilities), as shown in Fig. 10. Due to these 
instabilities, some fraction of the total energy is redi- 
rected into the kinetic energy of non-radial motion. This 
fraction is typically a few percent and can be high as up 
to 40 % locally in the transition layer between the hot 
bubble and the dense shell. Therefore less energy is avail- 
able for the kinetic energy of radial motion of the dense 
shell and the internal energy of the hot bubble. Another 
effect triggered by instabilities is the mixing of material 
in the hot bubble and the demise shell near the contact 
discontinuity which leads to a lowered temperature there. 
The consequent decrease in cooling time thus reduces the 
pressure and hence the extent of the hot bubble. 

The reduced extent of the hot bubble results in a re- 
duced luminosity £x,acis (by a factor of about 1.6, Fig. 
11). The luminosity in the medium energy bin (0.09- 
0.2 keV) remains almost unchanged. In the low energy 
bin (0.01-0.09 keV) the reduction of the extent of the 
hot bubble is compensated by the increased emission of 
lower temperature plasma resulting from turbulent mix- 
ing due to the instabilities, but we cannot quantify these 
effects accurately, because most of the energy in this en- 
ergy bin is emitted by only a few grid cells in our models 
(see §3.2.1). 
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Fig. 11. X-ray luminosity Lx.ACiS in the energy bins between 
0.09 - 0.2 keV and 0.2 - 10 keV as a function of the evolution time 
in the one-dimensional and two-dimensional runs B5 

5. COMPARISON WITH OBSERVATIONS OF APPARENTLY 
ROUND PNS 

In this section, we apply our one-dimensional mod- 
els to two PNs which show apparently round morphol- 
ogy. Rather than computing new 2D runs, we used ID 
runs in our modeling of BD-|-30°3639 and NGC 40, since 
other uncertainties, as, e.g., those due to uncertain abun- 
dances, are much larger than the differences between the 
ID and 2D runs. We use measured values of the fast 
wind properties (velocity, mass outflow rate) estimated 
from optical lines and derive the corresponding proper- 
ties of the slow wind (using eqns. 7 and 10) to fit the 
observed expansion velocity of the shell. In these simula- 
tions, we do not assume MfVf = MsVs- Then we choose 
the evolutionary age in our simulation, where the size 
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of the model dense shell is equal to the observed size of 
the PN, and compare the temperature of the hot bubble 
and the X-ray luminosity with values inferred from avail- 
able X-ray observations. Thus we constrain our models 
using four observed properties of PNs - t^exp, '''shcii, 
and Lx - whereas ASB06 take Tx, Lx and the dynami- 
cal age (rshou/wcxp) into account while comparing their 
models with the X-ray properties of PNs and ignore the 
measured values of the fast wind properties. 

5.1. BD+3(jP3639 

BD+30°3639 shows well-resolved, extended X-ray 
emission (Kastner et al. 2000) which lies within the in- 
terior of the shell of ionized gas seen in optical images. 
The extent of the shell in CHANDRA images is approx- 
imately 5" X 4" . Assuming a distance of 1 kpc (Kastner 
et al. 2000), the radius is therefore about 3.3 x lO'^^ cm. 
The fast wind velocity is about 700 km s~^ (Leuenhagen 
et al. 1996) and the mass loss rate of the fast wind is 
about 10"*^ M0 yr-i (Soker & Kastner 2003). The av- 
erage expansion velocity derived from [OIII] and [Nil] 
emission lines is 25.5 km . Kastner et al. (2000) give a 
CHANDRA ylCK spectrum of BD+30°3639 between 0.3 
- 1.7 keV and derive a total luminosity of 1.6 x 10'^^ org 
s~^ after fitting the spectrum with a variable-abundance 
MEKAL model. The deduced emission-region tempera- 
ture is 2.7 X 10^ K. 

A new simulation D7 with the observed parameters 
was performed using Ms = 7 x 10~^ yr~^. The 
observed size of the shell is reached at an age of 390 years. 
The model expansion velocity is 25.9 km s~^, which is in 
good agreement with the observations^. However, both 
the temperature of the hot bubble (1.1 x 10^ K) and the 
luminosity (9.2 x 10^^ erg s~^ in the full range of ACIS 
of 0.2-10 keV, 7.8 x 10^^ erg s~^ in the observed energy 
range) are too high in our model. 

The model spectrum (Fig. 12, top) is different from 
the observed spectrum. In the observations, there is a 
strong, broad peak at about 0.4 keV, which is not appar- 
ent in any of our model spectra. Furthermore there is a 
narrower, less prominent peak at about 0.9 keV in the 
observed spectrum, however, our model spectrum shows 
a strong, broad peak at about 0.9 keV (only visible in 
the models with a fast wind velocity above 600 km s^^). 
The strongest lines in the region around the peak at 0.9 
keV are iron lines (Fig. 12, top). This suggests that 
the iron abundance (and perhaps also those of other ele- 
ments) in BD+30°3639 are very different from the solar 
abundances (Anders & Grevesse 1989), which we have 
used in our models including run D7. 

5.1.1. The effects of uncertain abundances on X-ray 
emission properties 

The abundance problem may be especially important 
for objects where the central star is of late [WC] type, 
as in BD-|-30°3639 and NGC 40. In these objects, the 
abundances in the fast wind can be very different from 
those of the dense shell. The former, however, are very 
important for the X-ray properties, as the hot bubble 
consists of material ejected in the fast wind. As mixing 

It is interesting to note that Mf Vf = Ms Vg is satisfied in this 
object. 



of material from the dense shell into the hot bubble may 
occur, the abundances in the dense shell also may be 
relevant for the X-ray properties. Such mixing has been 
inferred in the hydrogen-deficient PN Abell 30 by Chu et 
al. (1997). 

An inspection of published values of abundances 
BD+30°3639 shows that the abundances, both in the 
fast wind and the dense shell, are in general poorly de- 
termined. Aller & Hyung (1995) estimated the following 
nebular abundances: C and N are solar, O is roughly 
solar, Ne is depleted (although they only give a lower 
limit) from a high spectral-resolution optical spectrum 
of BD-h30°3639. Bernard-Salas et al. (2003) determined 
that C is enhanced by a factor of 2, N and O are roughly 
solar, Ne is enhanced by a factor of 1.7, using IR spectra 
and UV observations. 

The abundances in the fast wind are determined by 
Leuenhagen et al. (1996) by fitting optical spectra. They 
find an upper limit of the hydrogen abundance and that 
helium, carbon and oxygen are very prominent as ex- 
pected. Arnaud et al. (1996) fit ASCA observations with 
two sets of abundances: (i) a solar C and half-solar N 
abundance, 0.2 solar for Ne, O and Fe are depleted; (ii) 
C is enhanced by a factor of 354, N and Ne are enhanced 
by a factor of 10, O is roughly solar and no iron. Kast- 
ner et al. (2000) concluded that C is enhanced, Ne is 
roughly solar and N and O are depleted, from a fit of 
its CHANDRA ACIS spectrum. Maness et al. (2003) 
revisited the CHANDRA observations of Kastner et al. 
(2000) and showed that the fast wind abundances (the 
second set of abundances in Arnaud et al. 1996) can 
fit the observations much better than the nebular abun- 
dances of Aller & Hyung (1995). Furthermore they find 
that an enhanced Ne abundance relative to solar is re- 
quired. Georgiev et al. (2006), re- investigating the same 
CHANDRA observations, also fit the spectrum with two 
sets of abundances: (i) C is enhanced by a factor of 3.7, 
Ne, N and O are depleted, only Si, Ca and Ni of the 7 
heavier elements are solar, the others (Mg, S, Ar, Fe) are 
not existent; (ii) C is enhanced by a factor of 350, N and 
Ne arc enhanced by a factor of 20, O is roughly solar, 
7 heavier elements up to Ni are solar. In summary, the 
available data can be fit with a wide range of abundances, 
although there seem to be general agreement that carbon 
is enhanced and oxygen is either depleted or solar. 

In order to demonstrate the dramatic effect of chang- 
ing elemental abundances on X-ray emission properties, 
we have recalculated the spectrum from our model D7 
using the elemental composition determined by Kastner 
et al. (2000) (i.e. only H, He, C and Ne, but with abun- 
dance ratios as in the Sun). The spectrum has the form 
as shown in Fig. 12 (bottom), which is closer to the ob- 
served shape. The total luminosity in the range of 0.3 
-1.7 keV is then reduced to 1.75 x 10'^^ erg s~^, which 
is now very close to the observed value of 1.6 x 10^^ erg 
s-i. 

However, an important caveat to note here is that the 
cooling function depends on the abundances. For exam- 
ple, at temperatures > 10® K the most important coolant 
is iron (see Fig. 18 in Sutherland & Dopita 1993). Hence 
an accurate calculation of X-ray properties as a function 
of abundances would require us to abandon the use of a 
cooling function, replace it with a set of rate equations for 
the important species and calculate the resulting cooling 
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Fig. 12. — X-ray spectrum for the one-dimensional runs D7 at an 
age of 390 years; Top: calculation with the default set of elements 
and identification of the strongest lines; bottom: here only the 
elements H, He, C and Ne have been taken into account 

for all transitions. This, however, would need a sub- 
stantially increase in the computational effort, which is 
beyond the scope of this paper. 

Such an effort, however, might be warranted for mod- 
eling the high resolution spectrum of BD-|-30°3639, re- 
cently taken using the LETG on CHANDRA, which 
shows abundance ratios of C/0 ~ 20, N/0 < 1, Ne/0 
~ 4 and Fe/0 < 0.1 (Kastner et al. 2006). The first step 
towards self-consistency of the cooling treatment would 
then be to modify the cooling function taking into ac- 
count the effects of lower metallicities as in Fig. 13 in 
Sutherland & Dopita (1993). 

5.2. NGC 40 

Another example of a PN, which is detected in X-ray 
and is almost round, is NGC 40 (Montez et al. 2005). 
Compared to BD-|-30°3639, the fast wind velocity is sig- 
nificantly higher in this object (1000 km s~^); the mass 
loss rate of the fast wind is 2.4 x 10~^ M0 yr^^ (Leuen- 
hagen et al. 1996). The observed radius of the shell is 
20,000 AU at an estimated distance of ~ 1 kpc. The av- 
erage expansion velocity from [OIII] and [Nil] emission 
lines is 27.5 km s~^. The measured X-ray luminosity is 
1.5 X lO'^" erg s"-'^ in the range of 0.3-1 keV and the in- 
ferred temperature of the X-ray emitting gas is 1.5 x 10® 
K. 

This object was modeled with a new run ElO. We chose 
a mass loss rate of the slow wind of 5 x 10~^ Mq yr~^ 
to adjust the expansion velocity of the dense shell in our 
model to the observed value^. The expansion velocity is 

* In this simulation, Mf V( is not equal to Ms f s , but lower by a 



27.03 km s^^, which is very close to the observed one, as 
expected by our choice of the parameters. The observed 
size of the bubble was reached in our model after 3775 
years, thus NGC 40 is much older than BD-|-30°3639. 

We find from our model that the temperature of the hot 
bubble is 2.2 x lO''' K and the luminosity is 2.4 x 10'^^ erg 
s-i in the full range of ACIS (0.2-10 keV) and 7.6 x 10^2 
erg s^^ in the observed energy range. If, as in the case 
of BD+30°3639, we restrict the composition to only H, 
He, C, N and O, the above luminosities are reduced to 
1.5 X 10"^^ erg s^^ and 6.4 x lO''^ erg s^^, respectively. 
These values are still too high (by a factor of 400) to 
explain the observations. 

5.3. The effects of different time histories of the fast 

wind 

The previous sections show that in our models of both 
objects, BD-|-30°3639 and NGC 40, and are higher 

than the observed values. Noting that is most sensi- 
tive to V[ (§3.1 and eqn. 11), we investigate whether bet- 
ter fits to the data could be obtained by assuming that 
the value of V[ was lower in the past (thus implying that 
both objects are now being observed in a phase where an 
even faster wind has emerged producing the line profiles 
analyzed by Leuenhagen et al. (1996) for deriving ff). 
We therefore compare two models, one with an ultra-fast 
wind of 1000 km s~^ (run BIO) and the second in which 
the fast wind has a velocity of 300 km s~^ for the first 
100 years followed by the ultra-fast wind (run B3-I-10). 
The mass loss rate of the slow wind was 7 x 10~® Mq 
yr~^. The momentum discharge MfVf = MgVs is again 
kept constant. 

We find that, in the second model, the contact discon- 
tinuity between the ultra- fast and the fast wind (CD2) 
reaches the original contact discontinuity between the 
fast and the slow wind (CDi) in 7 years and, after a 
short period of complex interaction, produces a hot bub- 
ble with a structure very similar to that in the other 
model (Fig. 13). The radial extent of the hot bubble is 
slightly smaller. The final expansion velocity of the shell 
is the same in both models (33 km s~^). 

We extract the luminosities at early ages in the evo- 
lution of these models (i.e., soon after emergence of the 
fast wind), when the size of the shell in each model has a 
value of 2000 AU. The corresponding ages are 340 (in run 
BIO) and 370 years (nm B3+10). The simulated lumi- 
nosities in the full range of ACIS are similar - 8.9 x 10^° 
erg s-i (run BIO) and 9.3 x 10^° erg s"^ (run B3-I-10). 
The temperatures of the hot bubble are 2.1 x 10^ K in 
both cases. Due to the slightly smaller extent of the 
hot bubble in run B3+10, the temperature gradient is 
slightly steeper. Although the spectra are very similar, 
a few small differences occur due to a larger quantity 
of cooler gas in the run B3+10: enhanced emission is 
seen in the following lines at 0.22, 0.29, 0.36, 0.43 and 
0.56 keV, which are fines of Si IX, C V, C VI, N VI and 
O VII. These results show that, once CD2 merges with 
CDi, neither the expansion velocity of the shell nor the 
total X-ray luminosity or spectrum distinguish between 
the very different time histories of the fast wind. 

Observable differences are only present in the short 
period of time between the emergence of the ultra-fast 

factor of 2. 
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Fig. 13. — Profiles of density (solid) and temperature (dotted) as a function of radius at ages of 100, 102, 104, 106 and 110 years for 
run B3+10 and for comparison at an age of 86 years for run BIO (bottom, right); the structure of the bubble is similar in both runs even 
directly after the merger 



wind and the merger of CD2 with CDi. We now have 
to check if it is likely that a specific PN is observed in 
this period. By using the radius of the dense shell at 
the time of emergence of the ultra-fast wind^ and noting 
that CD2 must expand faster than the fast wind, we can 
estimate a maximum time until the merger tcow- Thus 
for example, in our simulation B3+10, where the radius 
of CDi is about 6.3 x 10^^ cm after 100 years (Fig. 13), 



^coU 



G.3 X 10^ 



300 km s" 



6.65 years 



(12) 



We derive tcoii for BD+30°3639 and NGC 40, using the 
observed sizes of the two objects (3.3 x 10^^ cm for 
BD-|-30°3639 and 3 x 10^^ cm for NGC 40), and 300 

^ since the velocity of CDi is very small compared to that of 
CD2 



km s"'^ for the minimum velocity of CD2 since the tem- 
perature of the hot bubble is better fitted by such a fast 
wind velocity (cqn. 11). We find values of tcoU of about 
36 for BD+30°3639 and 320 years for NGC 40. Since 
both values are a small fraction (one tenth) of the total 
age of the bubble in each object, it is very unlikely that 
the observations of Leuenhagcn et al. (1996) caught these 
two objects in the special period after the emergence of 
the ultra-fast wind and before the merger of CD2 with 
CDi. 

In summary we conclude that the utilization of a 
smaller value of Uf in the past will not help in explain- 
ing the low values of and Tx observed in objects like 
BD+30°3639 and NGC 40. 

5.4. The importance of heat conduction and its 
inhibition by magnetic fields 
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Heat conduction across the contact discontinuity re- 
sults in a cooler and denser hot bubble. Zhckov & 
Perinotto (1996) give analytical estimates for the im- 
portance of heat conduction, calculating the ratio of the 
thermal conductivity time and the evolution time, claim- 
ing that heat conduction is important right from the be- 
ginning of the fast wind-slow wind interaction. However, 
they ignore the effects of radiative cooling, which are con- 
sidered by Soker (1994). By comparing the energy loss 
due to radiation and due to heat conduction from the hot 
bubble into the dense shell, Soker (1994) finds that radia- 
tive losses become significant compared to heat conduc- 
tion after a certain timescale tcooi (his eqn. 6), when the 
conduction front thickness (which increases with time) 
reaches a certain value. He derives this timescale by us- 
ing the growth rate of the conduction front from a model 
by Balbus (1986). However, both Soker (1994) and Bal- 
bus (1986) assume that the pressure of the hot gas is 
constant with time. Since the pressure of the hot bubble 
is expected to vary strongly with time in PNs (~ t~'^), 
the formulation of Soker (1994) may not be valid. 

The heat conduction front goes through three differ- 
ent stages. In the first stage, mass from the cold shell 
is evaporated into the hot bubble. This results in the 
hot bubble mass being dominated by the evaporating gas 
from the dense shell (Zhckov & Perinotto 1996). Thus 
in a PN in this stage the X-ray spectrum should reflect 
the abundances of the dense shell and not those of the 
shocked fast wind. This is also true for the second, quasi- 
static stage. In the third stage, i.e. "condensation", the 
region of intermediate temperature expands into the re- 
gion of the initially hot gas (Borkowski et al. 1990). In 
the last two stages radiative losses dominate heat conduc- 
tion. Soker & Kastner (2003) claim that in BD-|-30°3639 
the inferred depleted abundance of oxygen, based on the 
modeling of the CHANDRA ACTS spectrum by Kastner 
et al. (2000), suggests that the X-ray emitting conduction 
front is in one of the first two phases, i.e. located within 
the dense shell, making the implicit assumption that oxy- 
gen is depleted in the dense shell. However, we do not 
think that their claim is justified as the oxygen abun- 
dance seems to be roughly solar in the nebula, i.e. the 
dense shell (§5.1.1). The abundances as inferred from the 
new high-resolution X-ray spectra definitely show a very 
large overabundance of carbon (as expected in the fast 
wind from the [WC] central star of this object) in con- 
trast to the carbon abundance of the dense shell (§5.1.1). 
If heat conduction is operating in this object, it can not 
have evaporated a large amount of mass from the dense 
shell into the hot bubble. Alternatively, heat conduction 
is inhibited by magnetic fields in this object. 

Magnetic fields, which are believed to play an impor- 
tant role in coUimating the fast outflow, may reduce ther- 
mal heat conduction across field lines (e.g. Soker 1994). 
Following Borkowski et al. (1990), the conductivity coef- 
ficient can be written 

" ~ 1 + (n tan ^h/nh)^ 

U sing t he expression for the gyroradius S = 
\/2 Trrie/ (e B) with the electron mass nie and the mag- 
netic field strength B (Spitzer 1962), one can show that 
a very weak magnetic field of the order of 0.1 /iG is suffi- 
cient to reduce the conductivity (Soker 1994). Chevalier 



& Luo (1994) find in the models of magnetized stellar 
wind bubbles that the tangential magnetic field compo- 
nent in the shocked fast wind is larger than the radial 
component by several orders of magnitude, thus heat 
conduction is inhibited in radial direction. 

Heat flux inhibition by electromagnetic instabilities 
can also occur along the field lines, if /3 = Snp/B^ 1 
(Levinson & Eichler 1992), which gives a maximum mag- 
netic field 

B ^/8^Tp = VSnnkf = 




When the magnetic field is strong, the suppression is 

small (Pistinner & Eichler 1998). 

Several lines of evidence support the presence of mag- 
netic fields in PNs and therefore the inhibition of heat 
conduction is not unlikely. For example, magnetic fields 
in PPNs have been inferred from the presence of polar- 
ized OH and H2O mascrs (e.g. Vlemmings et al. 2006) 
and have been detected in the central stars of PNs (Jor- 
dan et al. 2005). Magnetic fields are also considered to be 
the main agent for producing bipolarity in PNs (Garcia- 
Segura et al. 2005). 

6. COMPARISON WITH ANALYTICAL RESULTS 

We compare our results now with the analytical re- 
sults of ASB06. A basic advantage of an analytical ap- 
proach is that it allows us to explore a large part of the 
parameter in a short time. However, the analytical ap- 
proach requires certain assumptions which at best may 
be only partially valid. Numerical simulations like those 
which we have carried out in this paper, although much 
more time consuming, are thus necessary for testing the 
limitations of the analytical approach as well as making 
realistic models of real objects. 

ASB06 calculated the structure of the hot bubble and 
the dense shell using self-similar solutions as in Chevalier 
& Imamura (1983). They computed the X-ray luminosity 
between 0.2 10 kcV (now called high energy bin) by tak- 
ing into account only those gas parcels in which the cool- 
ing timescale (from the cooling function of Sutherland & 
Dopita 1993) is larger than the evolution timescale of 
the flow. They derived the X-ray emissivities from the 
APEC database (Smith et al. 2001) and integrated the 
emission measure over the region between the inner shock 
and the contact discontinuity. Their results can now be 
compared with those from our simulations. 

We find a different temporal behavior of the luminos- 
ity in the high energy bin, Lx.AciSj compared to the re- 
sult of ASB06 whereas ix.Acis in our models increases 
slightly over the first few 100 years and decreases slowly 
(^ after that, in ASBOe's models it decays with 

time as ~ These differences probably arise because 
their cooling treatment is not fully self-consistent. They 
do not follow the evolution of the gas parcels with ra- 
diative cooling, as of course is done in our numerical 
simulations. Furthermore the temporal behavior of the 
luminosity in ASB06's models depends on the specific 
criteria for removing cooled regions (see their Fig. 9). 

The initial phase of low ix,Acis is of significantly 
smaller duration in our models as compared to ASB06 
- e.g. in model C3, it covers a range of only 100 years 
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in our models, whereas in ASB06 it is about 700 years. 
This difference is relevant to the understanding of X-ray 
emission from PPNs where the ages of the nebulae are 
typically ~ few x 100 years. The length of the phase of 
increasing luminosity is dependent on the parameters of 
our models the lower the velocity and the higher the 
mass outflow rate of the fast wind, the longer this phase. 
Both trends are qualitatively consistent with the results 
of ASB06. 

After the initial phase, a little step in the luminosity 
occurs (e.g. in nm B3 at an age of about 230 years in 
Fig. 3). Shortly before this transition, oscillations of 
density, pressure and velocity occur in the hot bubble. 
Frank & Mellema (1994) have also found oscillations in 
their numerical simulations and suggest that these occur 
when energies, e.g. the total and the internal energy, are 
almost equal. In our case, the internal energy is about 
96 % of the total energy. However, whether these os- 
cillations are only symptoms of, or the reason for, the 
transition, is unclear. Investigating the age at which the 
steps occur as a function of the parameters of the fast 
wind, we find that tstop ^ M^-^Wj ''^. Although we can 
not say as yet if the oscillations are a numerical arte- 
fact or a real physical effect, these scalings might help in 
eventually finding an answer. 

The temporal behavior of ivx,Acis hi the models of 
ASB06 and our models can be easily understood as fol- 
lows. In the self-similar models of ASB06, the radius of 
the contact discontinuity and that of the reverse shock 
are proportional to t, i.e. the volume of the hot bub- 
ble varies as t"^. The density decreases as t^^, hence the 
emission measure V varies a.st~^. As the temperature 
in the hot bubble remains roughly constant, the emissiv- 
ity A of the hot bubble gas is also constant and therefore 
the total luminosity (An^ V) varies as t~^. In compari- 
son, the empirical scaling laws in our runs are r^s ^ t'^'^^, 
red ~ V = (r3^ - ri) - t^.ss^ < p >^ ^-i-se and 

< p >2 y ~ resulting in Lx.ACis ~ The 

scaling laws of density and volume of the hot bubble im- 
ply that the total mass of shocked fast wind material 
increases as about ^ t^'^, whereas the mass input by the 
fast wind increases only as ~ t, since the mass outflow 
rate is kept constant. This difference in the time depen- 
dence reflects the evolution of the hot gas in the shocked 
fast wind under the influence of radiative cooling. In the 
beginning, only a small fraction of the mass injected by 
the fast wind is at high temperatures in the hot bub- 
ble, most of it lies near the contact discontinuity at low 
temperatures. This is because, due to radiative cooling, 
large parts of the shocked fast wind collapse towards the 
contact discontinuity. With increasing age, the density 
decreases and this effect becomes weaker, progressively 
increasing the fraction of gas injected by the fast wind 
which remains at high temperatures after being shocked. 

We find that the higher the mass outflow rate of the 
fast post-AGB wind, the higher Lx.ACiS is (§3.2), in ac- 
cord with the results of ASB06. However, within each 
set of our runs with constant M{V{, ix.ACiS decreases 
for increasing V{ - this behavior is different from the re- 
sults of ASB06 in which the models with V{ = 500 km 
s~^ (A5, B5, C5) always have the highest -Lx.acis value 
within each set of constant Mf vf. 



We can also compare the absolute luminosities in our 
models (which lie in the range 10^° - 10^"^ erg s~^) with 
those from ASB06® (sec Table 2). We find good agree- 
ment (i.e. within a factor of 2-3) for the runs in sets A 
and B, however, in set C the derived Xx,acis values of 
ASB06 arc higher (for example, a factor of 3 higher in 
run C5 and a factor of 5 higher in run C7). In addition, 
within each set (i.e. with constant MfVf), in ASB06 the 
values of -Lx.acis lie within a factor < 2 of each other, 
whereas in our models acis values cover a larger range 
of 4 10. 

ASB06 conclude that their models can fit the X-ray 
properties and the dynamical ages of specific PNs, in- 
cluding BD-f-30°3639 and NGC 40 (e.g. Fig. 5 in 
ASB06) . We think this conclusion is problematic for the 
following reasons: first ASB06 use smaller values of the 
fast wind velocities than those which have been derived 
from measurements (Leuenhagen et al. 1996). They find 
that then best fits to BD+30°3639 and NGC 40 arc ob- 
tained with models C3 and B4, respectively, for which 
the fast wind velocities arc too low compared with ob- 
servations (300 and 400 km s~^ compared to 700 and 
1000 km s~^, respectively), and justify their choice by 
claiming that the details of the late-time evolution of 
the fast wind are not important to the X-ray emission. 
They state that the X-ray emission predominantly em- 
anates from wind segments of the hot bubble which were 
expelled early, i.e. with a moderate velocity, and over 
a relatively short time (few x 100 years). Hence the 
X-ray emission of the hot bubble is determined by the 
moderate fast wind and the subsequent time history is 
unimportant; however, we find that the reverse is true 
- the X-ray emission is determined by the final stage of 
the time history during which the fast wind velocity has 
its largest value. 

Second, their claimed agreement between models and 
data appears to break down, when we consider in detail 
the model and observed values of all three parameters 
(Lx, Tx and age) together. Thus, ahhough BD-h30°3639 
is well fitted by 03 in the log(Lx)-age plane (Fig. 5 of 
ASB06), inthclog(Lx) log(Tx) plane (Fig. 6ofASB06), 
it is far removed from model C3 and lies close to model 
B6. NGC 40 (not shown in the log(Lx) age plane), is 
located at a position in the log(Lx)-log(rx) plane which 
appears to lie on an extrapolation of the curve for model 
B4 representing very late ages. However, taking into 
account the declining luminosity with increasing age and 
therefore size, this location would be reached at a size of 
the dense shell which is a factor 7 larger than observed. 

7. FITTING Tx AND Lx IN PNS - AN UNSOLVED 
PROBLEM? 

In section 5, we have found that using the values of 
t'exp, rsheii, Mi and Vf in BD-h30°3639 and NGC 40 as 
inferred from observations, our models produce values 

of Tx and Xx.ACiS which are much higher than the ob- 
served values. We explored three different mechanisms 
for reducing and ix.ACis- First, noting that the cen- 
tral stars in both objects are of [WC] type, deviations 

^ For this comparison, we have used bubbles with the same size 
rather than those with the same age, since the velocity of the con- 
tact discontinuity in our models is slightly different compared to 
ASB06. The velocities used by ASB06 follow from a cubic formula 
(Volk & Kwok 1985, also Table 1). 
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Fig. 14. — Plots of the mean density < p > of the hot bubble, the radius of the reverse shock Tts and the contact discontinuity r^d the 
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of the abundances from the solar vahies were tested as 
a way of reducing Lx,Acis- In the case of BD+30°3639, 
this mechanism appears to be promising, but needs to 
be tested using a cooling function consistent with the 
non-solar abundances. In the case of NGC 40, although 
ix,AciS was reduced, the discrepancy between the model 
and the observed value remains too large. The second 
mechanism, in which the fast wind had a slower velocity 
initially, was also not effective. Thirdly, heat conduc- 
tion, which has been discussed previously as a possible 
mechanism, docs lower Tx, but leads to a large increase 
in Lx.ACiS (Schonberner et al. 2006) due to an increase 
in the density of the hot bubble. We also found that, 
even though I/x,acis is reduced in models computed on 
a two-dimensional grid, the reduction is not adequate. 

Since Tx (Tx.acis) depends sensitively on V{ (Aff), is 
it possible that the values of these parameters inferred 
from observations by Leuenhagen et al. (1996) are too 
high? Since V{ is directly measured from Doppler shifts 
of absorption features, the uncertainties in determining 
Vf are likely to be small. The value of Mf is of course 
model-dependent and is also affected by the assumed dis- 
tance as Mf ~ d^/^ (Leuenhagen et al. 1996) - because of 
the distance dependency, Xx,ACiS should vary with dis- 
tance as ~ d^/^, since in our models ACIS ^ 
(Fig. 4). Lx.ACis in our model is also dependent on 
time and therefore on the size of the bubble. Assum- 
ing the dependency ix,AciS ~ t"'^ "'^ (see §6), we get 
i^XjACiS ~ r"" "^^ ~ Combining both dependen- 

cies, the ratio of Lx,ACiS in our models to the observed 
value varies as ^ cP-'^^~^-^'^ /tP = - because of 

this weak dependency, distance uncertainties cannot be 
invoked to reduce the discrepancies between the model 
and observed values of Lx.acis- Using the dependency 

-^x,Acis ^ Mf/"^ (Fig- 4), Mf would have to be reduced 
by a factor of 50 from the value inferred by Leuenhagen 
et al. (1996) in order to achieve consistency between the 
model and the observed value of Lx.acis in NGC 40. Of 
course, the problem of the model Tx being too high would 
still remain. In summary, the problem of hydrodynam- 
ical models for the formation of PNs producing values 
of Tx and Lx which are too high compared to observed 
values remains unsolved. 

8. CONCLUSION 

We computed the structure and evolution of, and X-ray 
emission from, numerical simulations of one-dimensional 
interacting winds models and compared the results with 
different analytical models. The kinematics and the evo- 
lution of our bubbles and the structure of the flows agrees 
well with analytical results. Comparing the X-ray prop- 
erties of the hot bubbles in our models with the ana- 
lytical results of ASB06, we find some agreements and 
many disagreements. The disagreements which are both 
qualitative and quantitative in nature argue for the ne- 
cessity of using numerical simulations for understanding 
the X-ray properties of PNs. 

We further investigated the luminosities in the energy 
ranges other than the 0.2 - 10 keV range covered by 
ACTS. We find that most of the X-ray flux may emerge 
at energies less than 0.09 keV (wavelengths above 137 A), 
although we do not have reliable quantitative estimates 
of this flux from our models. At these wavelengths, com- 



monly labeled as extreme and far UV, observations of 
PNs with FUSE (900 - 1200 A), EUVE (70 - 760 A), 
lUE (1150 - 3200 A) and GALEX for even longer wave- 
lengths (1350 2800 A) would therefore be quite impor- 
tant. Unfortunately, photospheric emission from the hot 
central star easily blends with the emission from the PN. 
Instruments with high spatial resolution are needed. 

We also investigated the X-ray spectra from our mod- 
els, which are not addressed by the analytical studies. 
We showed that the shapes of the spectra are diagnostic 
of the fast wind velocity which determines the tempera- 
ture of the hot bubble. 

We applied our spherical model to the objects 
BD-H30°3639 and NGC 40. We find in both cases that 
the simulated temperatures and luminosities of the hot 
bubble are much higher compared to the observations. 
While investigating these discrepancies, we identified sev- 
eral general issues which have to be considered for cor- 
rectly modeling the X-ray emission of PNs. 

The first one is the importance of an independent 
knowledge of the abundances - especially in objects with 
[WC]-type central stars, where the abundances of the 
hot bubble (shocked fast wind) can be dramatically dif- 
ferent from those in the dense shell (shocked slow wind). 
We could roughly reproduce the observed luminosity in 
BD-|-30°3639 by simply (i.e. in a non-self-consistent 
manner) removing several heavy elements including iron 
in the calculations of the X-ray properties. The treat- 
ment of cooling, however, has to be done self-consistently, 
as the cooling function itself is highly sensitive to abun- 
dances (e.g. iron). 

A second issue is the possible dependence of Tx and 
Tx on the time history of the fast wind velocity. Noting 
that Tx is most sensitive to Wf, we investigated whether 
Tx and Tx could be lowered by assuming that the value of 
V[ was lower in the past, but found that such a time his- 
tory will not help in explaining the observed low values 
of Lx and T^. The X-ray emission is determined by the 
final stage of the time history during which the fast wind 
velocity has its largest value. Although we have calcu- 
lated only one case of wind evolution with a particular 
jump in the fast wind properties, the claim of ASB06, 
that the X-ray emission of the hot bubble is determined 
by the moderate fast wind and that the subsequent time 
history is unimportant, is not supported by our results. 

We discussed the importance of heat conduction on 
lowering the temperature of the hot bubble, and its pos- 
sible inhibition by magnetic fields. As the presence of 
magnetic fields is supported both by observations and 
theoretical considerations, heat conduction is likely in- 
hibited. Even fields of the order of fiG are sufficient to 
inhibit heat conduction across (and along) the field lines. 
Since the X-ray spectra strongly indicate a very large 
overabundance of carbon in the hot bubble (as expected 
in the fast wind from the [WC] central star of this object, 
but much higher than the carbon abundance of the dense 
shell), if heat conduction is operating in this object, it 
can not have evaporated a large amount of mass from 
the dense shell into the hot bubble. Alternatively, heat 
conduction is inhibited by magnetic fields in this object. 

We also calculated a two-dimensional spherical model 
to examine the presence of instabilities and how they af- 
fect the structure of the hot bubble and the dense shell 
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and the resulting X-ray properties. Due to these insta- 
biUties, some fraction of the total energy is redirected 
into the kinetic energy of non-radial motion, which leads 
to slower expansion velocities of the shell (by about 20 
%) and a reduced luminosity between 0.2 - 10 keV (by 
~ 1.6), although the temperature of the hot bubble does 
not change significantly. 

Prom our detailed modeling, we conclude that the 
problem of hydrodynamical models for the formation of 
PNs producing values of and which are too high 
compared to observed values remains unsolved. The sim- 
ulations in our study are the first step towards modeling 
PNs and PPNs detected in X-rays which show bipolar 
or asymmetric structures with varying degrees of colli- 
mation. Our preliminary results show that higher colli- 
mation leads to lower X-ray luminosities (Stute & Sahai 
2006). Detailed results will be presented in a future pub- 
lication. 
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